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SUMMARY 


A previously developed local inviscid-viscous interaction technique 
for the analysis of airfoil transitional separation bubbles, ALESEP 
(Airfoil Leading Edge Separation), has been modified to utilize a more 
accurate windward finite difference procedure in the reversed flow region, 
and a natural transition/turbulence model has been incorporated for the 
prediction of transition within the separation bubble. Numerous 
calculations and experimental comparisons are presented to demonstrate the 
effects of the windward differencing scheme and the natural transition/ 
turbulence model. Grid sensitivity and convergence capabilities of this 
inviscid-viscous interaction technique are briefly addressed. A major 
contribution of this report is that with the use of windward differencing, 
a second, counter-rotating eddy has been found to exist in the wall layer 
of the primary separation bubble. 


INTRODUCTION 

The prediction of airfoil stall is an engineering problem of 
considerable importance as it determines one of the major operating limits 
of an aircraft. The physical mechanisms and different types of airfoil 
stall are reasonably well understood as evidenced by the excellent article 
by Tani (Ref. 1) which is nearly twenty years old. Nonetheless, the 
capability to predict stall has lagged considerably behind this 
progress due to the overall complexity of the flow, particularly the 
critical role played by the transition process from laminar to turbulent 
flow, and the development of reliable and accurate numerical methods for 
separated flow. It was pointed out by Tani (Ref. 1) that airfoils at 
moderate incidence angles, prior to either leading-edge stall or thin 
airfoil stall, experience local separation bubbles just downstream of the 
peak suction (minimum pressure) region. Therefore, as a start on the 
development of airfoil stall, the present investigation has been conducted 
to develop a method for the prediction of closed airfoil separation 
bubbles . 

Figure 1 shows a schematic diagram of an airfoil leading-edge bubble 
which occurs if the Reynolds number is sufficiently low so that the boun- 
dary layer remains laminar up to the minimum pressure point. Downstream 
of this point, separation occurs almost immediately since laminar boundary 
layers, in contrast with turbulent flows, are extremely sensitive to 
adverse pressure gradients • A separation bubble forms in which a 
recirculating streamline pattern is bounded by a shear layer. Since shear 
layer flows tend to be highly unstable to flow disturbances, transition 
from laminar to turbulent flow generally occurs in this shear layer. 

Further downstream, the turbulent mixing between the shear layer flow with 
the lower dead air region results in entrainment of higher energy air which 
energizes the flow near the surface thereby resulting in flow reattachment 
with subsequent turbulent boundary layer flow downstream. As shown in Fig. 

1 , the initial portion of the separation bubble is characterized by a 
pressure plateau followed by a pressure recovery region after the 



transition process is initiated, but prior to flow reattachment. An 
increase in incidence causes the bubble to move forward and contract in 
streamwise extent until the flow no longer reattaches. At this incidence 
angle, bubble bursting has occurred, thereby resulting in leading edge 
stall which is accompanied by an abrupt loss of lift since the suction peak 
has now collapsed with the resultant pressure distribution redistributed in 
a flattened form over the airfoil chord. 

There have been numerous experimental investigations conducted such as 
the work of Bursnall and Loftin (Ref. 2), Gault (Ref. 3), Gaster (Ref. 4), 
Horton (Ref. 5), Ntim (Ref. 6), Evans (Ref. 7), Roberts (Ref. 8), and more 
recently that of Mueller and Batill (Ref. 9) to provide information on the 
flow in the neighborhood of transitional separation bubbles. Owen and 
Klanfer (Ref. 10) deduced from experimental measurements that bubble 
bursting occurs if the momentum thickness Reynolds number at the laminar 
separation point is less than 125; similarly, Crabtree (Ref. 11) proposed 
a criterion based on the pressure rise over the bubble. Later, Horton 
(Ref. 12) developed a semi-empirical theory based on the experimental 
measurements of Gaster (Ref. 4) for the growth and bursting of laminar 
separation bubbles. At the present time airfoil analysis codes, such as 
the NASA-Lockheed multi-element airfoil code (Ref. 13) and the GRUMFOIL 
code (Ref. 14), use simple criterions such as these to detect the 
occurrence of laminar separation bubbles and whether or not bursting 
occurs. In these analyses, if laminar separation is detected without 
bursting then the flow is assumed to immediately undergo transition to 
turbulent flow with a jump discontinuity in the boundary layer parameters 
such as shape factor and skin friction. 

With the recent theoretical developments in boundary layer and 
inviscid-viscous interaction theory for separated flow, there have been 
several analytical investigations (Refs. 15-19) conducted to provide a more 
detailed description of the flow process in a laminar separation bubble. 
Further progress has been made by Gleyzes, Cousteix, and Bonnet (Ref. 20) 
in which their interaction analysis utilizes a natural transition model 
deduced from a correlation based upon amplification of laminar instability 
waves and free stream turbulence. 

The present report is the second in a series on the development of a 
prediction method for airfoil transitional separation bubbles. In the 
previous investigation, which was reported in Ref. 21 and subsequently 
condensed into the paper in Ref. 22, an inverse finite difference boundary 
layer procedure was iteratively combined with a Cauchy integral represen- 
tation of the inviscid flow which is assumed to be a locally linear pertur- 
bation to a known global viscous airfoil solution. Based on the favorable 
comparisons which were obtained with experimental data, the major conclu- 
sion drawn from this previous work was that this inviscid-viscous interac- 
tion model provides an accurate description of transitional separation 
bubbles provided the transition region is specified a priori. The focus of 
the current effort has been three-fold; first, to continue to improve this 
interaction technique through the inclusion of a windward differencing 
technique in the reversed flow region. Previously, the Reyhner and Flugge 
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Lotz (FLARE) (Ref. 23) approximation to the streamwise convection term was 
used in regions of reversed flow. Only minor changes have occurred in the 
pressure, skin friction and displacement thickness distributions with the 
inclusion of windward differencing. It has been found that a more accurate 
treatment of the convection term through the use of windward differencing 
has shown that a second, counter-rotating bubble exists in the interior of 
the main separation bubble. The existence of this secondary bubble 
structure appears to be insensitive to grid size as shown by numerical 
tests reported herein. In the second effort, the natural transition 
turbulence model of McDonald-Fish-Kreskovsky (Refs. 24 and 25) has been 
included to replace the forced transition model previously utilized. The 
McDonald-Fish-Kreskovsky transition/turbulence model has been found to 
predict transition too far downstream in separated and low free stream 
turbulent flows and hence, it is concluded from the present study that 
further development of this model is required to enhance its applicability 
to these flows. Finally, in a third effort, a brief study of the 
convergence of this interaction scheme has been performed along with 
additional comparisons with experimental data to further evaluate the 
scheme. Numerical results indicate that the present inviscid-viscous 
interaction technique is capable of reducing residuals to desired levels. 

LIST OF SYMBOLS 


a Structural coefficient 

c Airfoil chord 

Cf Skin friction coefficient 

Damping factor applied to mixing and dissipation lengths 
f Perturbation stream function 

F Velocity ratio, u/u g 

g Total enthalpy ratio, H/H g 

H Total enthalpy 

H Mixing length or ratio of local to edge density * molecular 

viscosity product 

L Reference length or dissipation length 

m Perturbation mass flow 

n Coordinate normal to reference displacement surface 

N Coordinate measured normal to reference displacement surface from 

the body surface 
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Pr 


Prandtl number 


\S 


Pr^ Turbulent Prandtl number 

q Magnitude of fluctuating velocities 

Re Reference Reynolds number 

Reg Local momentum thickness Reynolds number 

Rq Correlated momentum thickness Reynolds number 

R t Turbulent Reynolds number 

s,S Coordinates along reference displacement surface 

Tu Turbulence level 

u Velocity component parallel to reference displacement surface 

v Velocity component normal to reference displacement surface 

V Transformed normal velocity in Prandtl transposition theorem 

a Windward differencing weighting operator 

3 Pressure gradient parameter 

<5 Boundary layer thickness 

6* Displacement thickness 

6 t Stress thickness 

£ Eddy viscosity coefficient 

k von Karman constant 

D Transformed normal coordinate 

V Kinematic viscosity coefficient 

y Molecular viscosity coefficient 

£ Transformed tangential coordinate 

4> Velocity potential 

P Density 

^ Stream function 
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U) 


Interaction relaxation parameter 


Subscripts 

e Edge of boundary layer 

I Inviscid 

ref Reference solution 

tj Start of transition 

End of transition 

T Turbulent 

v Viscous 

00 Free stream 

1 Start of interaction region 

2 End of interaction region 
Superscripts 

1 Perturbation quantity 

+ Inner wall non-dimensionalized coordinate 

k Global inviscid-viscous iteration counter 


INVISCID-VISCOUS INTERACTION ANALYSIS 

Since the initial development (Refs* 21 and 22) of the transitional 
inviscid-viscous interaction technique, work has continued to improve this 

analysis for short transitional snpaxatjLoi Lbubbl ^s * In this analysis, the 

two-dimensional boundary layer equations are solved in inverse form (Ref* 
26) iteratively with an incompressible Cauchy integral perturbation 
analysis for the inviscid flow* Interaction between the inviscid and 
viscous solutions is accommodated by using an update formula (Ref. 27) 
which modifies the specified displacement thickness based upon the 
differences between the inviscid and viscous velocities at the edge of the 
boundary layer. 

In this section, a review of the inviscid analysis, boundary layer 
analysis, and the interaction iteration procedure are presented along with 
a detailed description of the windward differencing scheme and the natural 
transition/turbulence model which have been used in the present study. 
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Inviscid Analysis 


The local inviscid analysis which was used in the present investiga- 
tion as well as our previous work (Refs. 21 and 22) assumes that the 
disturbance field induced by the presence of a transitional separation 
bubble can be treated as a small disturbance to the global airfoil flow. 

In Appendix A, a perturbation analysis is used to deduce the governing 
equations and boundary conditions which describe the local disturbance to 
the inviscid flow due to the presence of a protrusion on an airfoil. In 
this analysis , the inviscid flow is split into an inner region which is in 
the immediate vicinity of the protrusion and an outer region which 
describes the global flow about the entire airfoil. Under a particular 
limiting condition for the protrusion height to length ratio, it is shown 
that to lead order the inner inviscid flow problem reduces to that of a 
protrusion placed on a flat surface subject to a uniform incoming stream. 
The magnitude of this uniform stream velocity is the airfoil local surface 
speed deduced from the global or outer airfoil solution without the 
presence of the protrusion. It is shown in Appendix A that the governing 
equation for the perturbation potential for the inner flow region is 
Laplace's equation subject to the usual small disturbance surface boundary 
condition at the airfoil surface and uniform flow at large distance from 
the bump. It is concluded on the basis of this analysis that to lowest 
order the disturbance field created by the displacement thickness induced 
by a transitional bubble can be represented by a lineal source distribu- 
tion placed on the airfoil surface at the transition site. In the future, 
this technique could be extended to include compressibility effects and 
thereby develop a rational theory for the disturbance field induced by an 
airfoil separation bubble in compressible flow. 


In Ref. 21 a detailed discussion was presented which showed that if a 
lineal source distribution is used to represent the local inviscid flow 
then the surface speed is given by 


u'(s,o) 




(i) 


The source strength is proportional to the streamwise rate of change of the 
difference between the computed displacement surface and the reference dis- 
placement surface both of which are shown schematically in Fig. 2. 

The Cauchy integral given in Eq. (1) is evaluated from s, to s„ which 
is the region of strong interaction that results from the presence of the 
local transitional separation bubble. In the present calculations, the 
reference displacement surface is computed from the GRUMFOIL analysis (Ref. 
14) in which instantaneous transition from laminar to turbulent flow is 
assumed to occur at the predicted laminar separation point. As a result, 
the interaction displacement surface merges smoothly with the reference 
displacement surface upstream of the interaction as shown in Fig. 2. 

Hence, the lower limit, s^, on the Cauchy integral is placed sufficiently 
far upstream of the interaction region where the source strength is zero. 
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The placement of the downstream limit is somewhat more complicated since 
the interacted displacement surface does not merge into the reference dis- 
placement surface downstream of the interaction region. However, this mis- 
match does not present a problem since, as was demonstrated in Ref. (21), 
the interacted and reference displacement thicknesses typically are nearly 
parallel to each other and therefore, do not contribute to the Cauchy 
integral in this downstream region. 

The Cauchy integral given in Eq. (1) is evaluated numerically using a 
second order scheme developed by Napolitano, et al. (Ref. 28). This 
numerical procedure permits the use of a nonuniform mesh distribution which 
was used in the present problem to concentrate grid points near the center 
of the interaction region in order to adequately resolve the high gradient 
phenomena which occurs as the flow undergoes transition from laminar to 
turbulent flow. The same surface grid point distribution was used in both 
the inviscid and boundary layer analyses thereby avoiding interpolation 
between these two solutions in the interaction calculations. 

Viscous Analysis 

The viscous solution technique used in the present investigation is 
the inverse boundary layer procedure presented by Carter (Ref. 26) for the 
analysis of separated flows. Although the inviscid analysis discussed in 
the previous section has been limited to low speed flow, the boundary layer 
analysis used in the present study was adapted from earlier work which was 
for compressible flow. Thus the fully compressible boundary layer analysis 
is presented here. 

Formulation 


The nondimensional boundary layer equations are written as follows in 
terms of the reference displacement surface coordinate system shown 
schematically in Fig. 2: 



The v-comp onent of velocity and the n-coordinate are scaled in the usual 
manner by /Re^c where Reoo C is the Reynolds number based on the free 
stream flow conditions and the airfoil chord. The boundary conditions 
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imposed on these governing equations are: 


dH 

n = -S* ef u-v=0, h or Specified 


( 5 ) 


n— oo u-— u e i h — H e 

The Reynolds stresses are related to the mean flow by 


du 


-P uV , -p v'h' 


< y dn 


Pr t dn 


( 6 ) 


(7) 


where £ is the viscosity coefficient. The transition from laminar to 
turbulent flow may be modeled using a streamwise intermittency factor, 

Y(s), which varies from 0 to 1 over a specified region. 

Werle and Verdon (Ref. 29) showed that it is convenient to transform 
the boundary layer equations with the Prandtl transposition theorem which 
is given by 

* d &ref 

S =S , N^n + Sref « V J v + u ^ (g) 

where N is a transformed normal coordinate measured from the airfoil 
surface perpendicular to the reference displacement surface. In Ref. 29 it 
is shown that the form of the boundary layer equations is unchanged by this 
transformation and the same boundary conditions (Eq. (5)) are imposed at 
N = 0. 


The development of the inverse formulation begins by transforming the 
equations, expressed in primitive variables, by the following transforma- 
tion of the independent variables 

N 

i -f P 9 ^e u e ds 17 = S 7 / dN 

0 0 

which is quite similar to the Levy-Lees transformation. It is helpful to 
scale the normal coordinate by the displacement thickness in strongly 
interacting flows since this step insures that the boundary layer thickness 
is approximately constant in the transformed coordinate. The continuity 
equation is eliminated by introducing the stream function 


P u = 


d± 

dN 


P* = " 


d± 

dS 


( 10 ) 
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The value of the stream function at the boundary layer edge is written in 
terms of the displacement thickness 

♦ — />e u e (N “ (11) 

where 


8 *=/° 

*0 ^e u e 


Then with the definitions 


m = 


= Ve 8 * h = / (~f ” ') dT ? 


the edge value of the stream function can be written as 


+ m (tj- | + h) os r, — co 

A perturbation stream function is defined as 


( 12 ) 


(13) 


(14) 


f = ~ (*? ~ • + h)| with F s (15) 

such that f 0 as n -*■ 00 for a prescribed m. Note that in the transformed 
inverse formulation the perturbation mass flow, m, is prescribed and not 
just the displacement thickness. 


Transformation of the compressible boundary layer equations with Eq . 
(9) and the introduction of the perturbation stream function defined in Eq. 
(15) results in the following set of governing equations: 


<?T _ m 
d V ” sfzl 


( I - 7) - h) 4 1 

drj 


(16) 


m‘ 


F ♦ mF to + h j|^=m 2 /9(g-F 2 ) 


(17) 


m 2 F - m [,/z£ T + mF - | + h)] 


Pr to) L \ ^ Pff/ dr, j 


( y J±2L- A r f /, I\ F 


ar 

dij 


dr, 


(18) 
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where 


g 



_|_ dMe t PH- 

^"*ed i 'Wt 


(19) 


Equations (16)— (18) are solved for F, g, f and E for a prescribed 
streamwise distribution of m subject to the following boundary conditions: 


i) - 0 


77 — CD 


F=r=o 9 =g„ or 

\ 

F = g — ► 1 and f 


specified 


( 20 ) 


These equations can also be solved in the direct mode with 6 prescribed and 
the outer boundary condition f = 0 eliminated. In this case, if m is set 
equal to /25 then the usual Levy-Lees formulation is deduced with the only 
difference being that the normal component of velocity has been 
re-expressed in terms of the stream function. In the inverse case, the 
unknown pressure gradient parameter is deduced simultaneously with the 
remainder of the solution. The numerical solution of these equations for 
the direct and inverse mode is an implicit finite difference technique 
which is first order accurate in the stream direction and second order 
accurate in the normal direction. The details of the numerical scheme are 
presented in Ref. 26. In the next section, the modifications made to this 
scheme in the reversed flow region are discussed. 

Windward Differencing 

In the previous formulation (Ref. 21) of the finite difference repre- 
sentation of the boundary layer equations, the FLARE (Ref. 23) approxima- 
tion was used in regions of reverse flow (i.e., within the separation 
bubble) to insure numerical stability. This approximation assumes that 
streamwise convection is zero in these regions. An objective of the cur- 
rent study has been to assess the accuracy of the FLARE approximation with 
a more accurate "windward" finite difference scheme in which the flow 
direction is properly accounted for in the finite difference representation 
of the streamwise convection terms. 


The finite difference cells shown in Fig. 3 are used in the numerical 
approximation to the boundary layer equations given in Eqs. (16)— (18) for 
attached flow. A second order accurate central difference operator is 
applied to the normal gradient terms whereas a one-sided first order 
accurate operator is applied to the streamwise gradient terms. This first 
order accurate streamwise differencing scheme is known to be much more 
stable and free of oscillations than the second order Crank Nicolson 
difference scheme due to the numerical damping (diffusion) inherent in the 
lower order accurate scheme. 
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Since the boundary layer equations are parabolic, an instability will 
arise when the solution inarching direction is opposite to the flow 
direction. Reyhner and Flugge Lotz (Ref. 23) have shown that this 
instability is easily avoided by assuming that the streamwise convection 
terms are zero in reversed flow regions. For the present equations, the 
FLARE approximation becomes 


when F < 0. 


^ ^ = 0 
dt dZ 


( 21 ) 


This approximation is particularly attractive since special 
differencing cells do not have to be utilized in the implementation of the 
FLARE approximation. Therefore, the differencing cells shown in Fig. 3 may 
be used everywhere within the boundary layer thus simplifying the numerical 
scheme. 

It is apparent, however, that when this approximation is used, the 
converged solution of the numerical difference equations is not necessarily 
the solution to the boundary layer equations since a loss of accuracy due 
to the negligence of the streamwise convection terms is incurred. As an 
improvement to the FLARE approximation, a windward finite difference cell 
such as that shown in Fig. 4 can be used to calculate streamwise gradient 
terms in reversed flow regions. The streamwise gradients for both attached 
and separated flow are represented by the weighted difference operator: 


dF 



a 




+ 


( i -a) 


k- I 

F i+i,i 




( 22 ) 


where a = 1 F ,> 0 attached flow 
= 0 F < 0 reversed flow 
k = global iteration index 


The weighting coefficient, <*, is usually 1.0 which yields the usual first 
order backward difference operator for attached flow. For negative 
streamwise velocities, a is changed to 0 and the difference operator 
switches to a first order forward difference. Since the marching direction 
of the numerical scheme is always in the positive (i.e., increasing i) 
^-direction, the information used to calculate the forward difference 
operator of Eq. (22) must come from the previous global iteration results. 

The operator given by Eq. (22) may be applied to all of the streamwise 
gradient terms including the stream function ^-gradient in Eqs. (17) and 
(18). Numerical calculations, however, have indicated that the use of this 
operator on the stream function streamwise gradient leads to an instability 
as the bubble size increases or as the streamwise grid spacing becomes 


11 



J 


smaller. This instability appears to originate in regions where the normal 
velocity component is high such as exists at the ends of a large 
recirculating region. Application of the weighted difference operator on 
the stream function has been found to amplify the normal velocity component 
in these regions. Similar sensitivity to windward differencing on this 
term was observed in Ref. 30 in which a stream function— vorticity formu- 
lation of the boundary layer equations were solved for separated flow. 
Further research is needed to determine a stable iterative strategy if the 
d\p/d£ term is to be treated in the same manner as the 9F/35 term. Hence, 
this study has focused on the application of the windward difference 
operator to the streamwise convection term only. 



Natural 


Transition/Turbulence Model 


In the previous interaction work presented by Carter and Vatsa (Refs. 
21 and 22), the Cebeci-Smith turbulence model (Ref. 31) was used to deter- 
mine the turbulent eddy viscosity coefficient, e. Either an instantaneous 
transition or the intermittency function empirically developed by Dhawan 
and Narasimha (Ref. 32) was used to force the flow from laminar to 
turbulent. In this earlier report, with the proper selection of transition 
point and transition length, the effects of the separation bubble on the 
inviscid loading distribution could be predicted with the interaction 
technique. Unfortunately, the values of these quantities are not usually 
known a priori and the solution is sensitive to their specifications. 


A focus of the current effort, therefore, has been to incorporate a 
turbulence model into the interaction technique which has the ability to 
predict transition automatically. Briley and McDonald (Ref. 15) reported 
success in predicting transitional separation bubbles with their time 
dependent Navier-Stokes/boundary layer solver using the McDonald-Fish 
turbulence model (Ref . 24) . This model is based upon the solution of the 
integral form of the turbulent kinetic energy equation. Since this model 
has been presented a number of times in previous reports, including the 
recent critique given by Walker and Werle (Ref. 33), only a brief discussion 
sion will be presented here. The turbulent kinetic energy equation is 


where 



(23) 


(24) 



d u/u® 3 

t 1 


1 s_ 

drj 

— 6t) , 

(25) 
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( 26 ) 



d _ u/Ue \ + 

dr, ) 





E = 


dS r 

^e U e ds v c 


where 


f T = I - COS 


7T N 

§7 


(27) 


(28) 



McDonald - Fish Model 
McDonald -Krtfkovtky- Fish Model 


(29) 


In the numerical predictions to be presented in the Results Section, the 
McDonald-Fish-Kreskovsky model (Ref. 25) has been used in attached direct 
boundary layer calculations and the McDonald-Fish (Ref. 24) model has been 
used in separated interacting calculations as suggested by Briley and 
McDonald. For the McDonald-Fish-Kreskovsky model, the value of 6 T is 
defined as the first N location from the edge of the boundary layer where 
the local shear stress exceeds 2 percent of its maximum value. If 6 X is 
computed to be less than the boundary layer thickness, 6, then 6 T is reset 
to equal 6. For the McDonald-Fish model, <5 T is set equal to 6. Structural 
coefficients, a^ , a 2 S and a^ , have been introduced in terms of the local 
mixing length, £, to determine the components of the perturbation 
velocities . 


u' v' 


= a. 




— 2 
T 


}■/ 


d u 

d N 


d u 
<?N 


(30) 


u 


/ 2 


= Q2 Q 2 


(31) 



The total magnitude of the fluctuating velocities is 


(32) 
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q 2 = u' + v' + w' 


(33) 


The local mixing length model as given by McDonald and Fish (Ref. 24) and 
later modified by McDonald and Kreskovsky (Ref. 25) is 


A - % 1& 

8 ' S 



tanh 



for y < S 



y > S 


(34) 


Similarly, the dissipation length is 


kr =2> -?• tanh 



(35) 


where k = 0.43 is the value used for the von Karman constant and is the 
damping factor given by 




+ erf 


f N + 

^ O l 


— }] 


By/~2 


The inner region scaled normal coordinate is 


(36) 


N + = N y r/p 

v 

The free stream dissipation length is 



(37) 


(38) 


The remaining unknown quantities are the values given to the structural co- 
efficients, a^ , a 2 » and a^ • McDonald and Fish showed that a, is nearly 
zero in laminar flow and reaches a value of 0.15 for fully turbulent 
boundary layers. McDonald and Fish derived a relationship between a-^ and 
the turbulent Reynolds number, R T , where: 


J. f^T 

S T J 0 >'T dN 

T ” , -Sj (39) 

JjJo vdN 


where v and v^, are the laminar kinematic and turbulent eddy viscosities, 
respectively, and 6^. is inner layer thickness defined as the first location 
from the wall where v is approximately 4 percent of the total effective 
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viscosity, v + v~. By first converting R T into a correlated momentum 
thickness Reynolds number, Rq, using an empirical relationship 

Rg « 100 R t 0 22 for R T < I , 


= 68.26 R T - 614.33 forR T >40 


( 40 ) 


the a^ structural coefficient may be calculated from the relationship: 


Q 


where a = 0.0115 
Rq = 100 




1+6.666 a 0 



(41) 


O 

A cubic polynomial which matches values and derivatives at both ends is 
used to express R^ in terms of R t for 1 <_ R T < 40. 


Briley and McDonald (Ref. 15) reported that their solution in separa- 
ted flows was sensitive to the magnitude of the normal stress terms in Eq. 
(23). Without a modification to structural coefficients, a£ and a^ , the 
flow upon separating, would not reattach and their solution diverged. For 

this reason, they allowed the difference in an< ^ a 3 appearing in Eq. (26) 
to vary linearly with R T using 


°2 “ °3 ” 0.3 + 0.6 (l - R r ) for R r < I 
= 0.3 for R r > | 


(42) 


Upon substitution of Eqs. (24)-(42) into the turbulent kinetic energy 
equation, the free stream mixing length, £«>, may be calculated. The local 
turbulent eddy viscosity is finally calculated from 


= 


/ / 
p U V 

~7u 

<3N 


pi 2 ^ 
d N 


(43) 


The transition between laminar and turbulent flow is primarily controlled 
by the source terms in Eq. (23) and by Eqs. (39)-(41) which determine the 
magnitude of the structural coefficient a^ . The source term E which 
controls the growth of the turbulent kinetic energy is proportional to the 
free stream turbulence level given by q^. For convenience q^ is normalized 
by the boundary layer edge velocity u^ as 


T 



I 


u e 


(44) 
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Interaction Iteration Procedure 


The present analysis is based on a inviscid-viscous iteration tech- 
nique which was previously presented by Carter (Ref. 27) and is adapted to 
the present investigation as outlined in Fig. 5. This procedure, which has 
been referred to as a semi- inverse technique by LeBalleur (Ref. 34), 
combines an inverse boundary layer technique with a direct inviscid analy- 
sis via the update procedure shown in Fig. 5. The key feature of this 
iteration procedure is the simple update formula 



which permits an inverse boundary layer analysis to be directly linked to 
an inviscid analysis which accounts for displacement thickness effects. 
Note that for simplicity the update formula is shown in Fig. 5 with the 
relaxation factor, oj > set equal to one. It was found by Kwon and Pletcher 
(Ref. 17) that convergence could be accelerated by making several inner- 
loop passes through the Cauchy integral and the update formula with the 
boundary layer prediction of the edge velocity, u g , frozen at its current 

global iteration value. This technique was used in the present calcula- 
tions with three inner-loop passes and was found to accelerate the global 
convergence rate by a factor of three with a 50 percent reduction in 
computer time as compared to calculations made without this inner 
iteration. 


RESULTS 

Windward Differencing 

In this section, the results obtained with the ALESEP (Airfoil Leading 
Edge Separation) code using the windward differencing operator of Eq. (22) 
in the reversed flow regions are compared with the results presented by 
Carter and Vatsa (Ref. 21) using the same code but with the FLARE 
approximation. In all of these calculations, a total of 100 grid points 
were used in the normal direction in the boundary layer analysis with these 
points distributed so that the minimum spacing was at the wall. The 
Cebeci-Smith turbulence model was used with the start and length of 
transition imposed at the same values as those reported by Carter and Vatsa 
with the use of the Dhawan and Narasimha intermittency function. 
Calculations using the windward differencing and FLARE approximations are 
shown along with experimental data for the Gault (Ref. 3) NACA-0010 airfoil 
and the Gaster (Ref. 4) series I-IV experiments. 

Gault NACA-0010 Airfoil 


The particular Gault case (Ref. 3) which has been analyzed here is for 
an NACA-0010 airfoil at an 8 degree angle of attack and a chord Reynolds 
number of 2 x 10 6 for which a transitional separation bubble occurred in 


16 





the leading edge region. This case has been chosen for the present 
assessment since the maximum reversed flow velocity, u/u = -.28, was the 
largest of any case analyzed in the previous investigation (Ref. 21). 

Results obtained from the GRUMFOIL code (Ref. 14) have been used for the 
reference surface velocity solution and reference displacement thickness 
distribution. A total of 71 grid points were distributed in the streamwise 
direction extending between s/c = 0 to s/c = .32 with the minimum spacing 
located in the transition region. The onset of transition was located at 
s/c = .0283 with a transition length of .0161. 

Figure 6 shows the predicted distributions of pressure, skin friction, 
and displacement thickness using the windward and FLARE approximations. 
Comparison of the pressure distributions shows that the only noticeable 
difference between the windward and FLARE predicted results is a slight 
decrease in the pressure coefficient at the "breakpoint" (transition point) 
with the windward scheme. A similar expansion at the end of the pressure 
plateau region has been observed in the recent experimental data of Jansen 
and Mueller (Ref. 35) as well as the earlier work by Horton (Ref. 5) on 
finite swept transitional bubbles. Comparison of the computed results with 
the experimental pressure data in Fig. 6(a) shows that the inability of the 
analysis to predict the constant pressure region near the peak suction 
pressure is still unexplained as it is clearly not affected by the improved 
differencing procedure used in the reversed flow region. Figures 6(b) and 
6(c) show that, in general, only small differences exist between the 
computed skin friction and displacement thickness distributions due to the 
inclusion of the more accurate windward flow differencing technique. Hence, 
it is concluded from the analysis of this case, which contains a 
significant backflow velocity of u/u g = -.28, that the FLARE approximation 
is quite accurate in predicting the overall effects of the transitional 
separation bubble. Comparison of the windward and FLARE results shows, in 
this case, only small differences which principally occur in the reversed 
flow region. These differences, though, are very interesting as the 
inclusion of windward differencing has revealed a new structure in the 
recirculating flow region. This change in bubble structure is discussed 
below. 

It is observed in the skin friction distribution in Fig. 6(b) that the 
calculation performed with the windward scheme in contrast with the FLARE 
result shows a small region of forward flow (C^ > 0) occurs in the 
interior of the separated flow region. This region occurs at the same 
location where the surface pressure predicted with the windward scheme 
shows a slight expansion at the end of the plateau region. Figure 7 shows 
a comparison of the computed streamlines in the viscous region obtained 
with the windward and FLARE schemes. Overall, these streamline patterns 
are very close with the major difference being that the more accurate 
treatment of the reversed flow region, via the windward scheme, has 
revealed the existence of a second, counter-rotating bubble inside of the 
primary separation bubble. To our knowledge this is the first time such 
a structure has emerged from numerical calculations of the interacting 
boundary layer equations for closed separation zones on solid surfaces. 
Physically, such a structure is known to exist in separated flows as 
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evidenced by several figures in the excellent compilation on flow 
visualization recently published by Van Dyke (Ref. 36). 

Velocity profiles at eight locations for the windward and FLARE calcu- 
lations are shown in Fig. 8. with their positions in and around the bubble 
region denoted by the arrows in Fig. 7. The profile at s/c = .031 shows 
the small region of forward flow near the wall which is inside of the 
secondary bubble predicted by the windward scheme. Except for the region 
where this inner bubble exists, use of the windward differencing scheme 
results in slightly lower local velocities within the boundary layer than 
those calculated by the FLARE approximation. This difference accounts for 
the slight increase in the displacement thickness found in the windward 
scheme over that deduced by the FLARE technique for this case. 

A grid study has been conducted to determine the double bubble 
structure sensitivity to the streamwise grid spacing. Only the streamwise 
grid was varied since the normal mesh of 100 grid points across the boun- 
dary layer was thought to be adequate to resolve this secondary bubble. 
Within the same calculation region (0 < s/c < .32), a 31 point and a 141 
point grid were used in the windward and FLARE schemes. As for the 71 
point grid discussed above, the points for these two calculations were 
distributed nonuniformly with the minimum spacing in the transition region. 
Figure 9 shows the effect of streamwise grid spacing on the skin friction 
coefficient using the windward scheme. Despite some differences in the 
skin friction due to the change in the mesh size, the presence of the inner 
bubble remains unchanged for these calculations thus demonstrating that its 
existence is not grid sensitive. Similarly, it was found that the 
structure of the bubble using the FLARE approximation did not change from 
that shown in Fig. 7 with these modifications to the streamwise mesh. 

Convergence histories for the NACA-0010 calculations using the wind- 
ward and FLARE schemes are shown in Fig. 10 for the three streamwise grids. 
A relaxation factor of 0.5 was used in Eq. (45) for both the windward and 
FLARE calculations for the 141 point grid calculation in order to overcome 
an iteration instability. A similar situation arose in the earlier work of 
Carter and Vatsa in which it was found necessary to use underrelaxation on 
the perturbation mass flow parameter, m, in the update procedure to 
eliminate a similarity instability. No attempt was made in the 141 grid 
point calculations to optimize the relaxation factor to obtain the fastest 
convergence rate. Relaxation factors of 1.0 were used in the 31 and 71 
point grid calculations. Figure 10 demonstrates that the convergence rate 
of the windward scheme is nearly the same as that for the FLARE and that 
both schemes converge monotonically to the desired level of residual. 
Windward differencing can, therefore, be used with little effect on the 
convergence rate of the iteration scheme if the fine details of the 
reversed flow region are of interest. Experience^has indicated that 
convergence of the edge velocity residuals to 10 produces sufficiently 
accurate results. 
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Gaster Experiment 


A second case with a less intense backflow region has been investiga- 
ted to compare the results obtained with the windward and FLARE 
differencing, and in particular, to see if the secondary bubble occurs with 
the inclusion of the more accurate windward differencing scheme. An 
experimental investigation of a transitional separation bubble on a flat 
plate performed by Gaster (Ref. 4) was chosen for this case. The 
separation bubble was induced by the pressure field generated by the 
placement of an inverted airfoil near the plate. Interaction calculations 
using the windward and FLARE approximations were performed for the series 
I- IV experiment of Gaster. The reference pressure distribution which was 
used in these interaction calculations was that measured by Gaster on the 
flate plate with the boundary layer tripped to turbulent flow near the 
leading edge. The reference displacement thickness was found from a direct 
boundary layer calculation with the reference pressure distribution 
prescribed and transition forced to instantaneously occur at the experi- 
mental trip location. A uniform mesh consisting of 81 streamwise grid 
points between s/L == 0.5 and s/L = 1.5 was used in the calculations. The 
Cebeci- Smith turbulence model was used with an assumed instantaneous 
transition at s/L = 1.0. 

Figure 11 shows only minor differences between the windward and FLARE 
difference schemes for the skin friction and displacement thickness dis- 
tributions for this case. There are no visible differences in the predicted 
pressure distributions of the two calculations. Also, the structure of the 
separation bubble does not change from that predicted with the FLARE 
approximation in this case. Figure 12 shows a comparison of the velocity 
profile for the two calculations located at s/L = 1.0125 where the magni- 
tude of the reversed flow velocity is a maximum. Only minor differences 
exist between the two profiles for this case. The local velocity reaches a 
u/u e = -.15 in the reversed flow region whereas in contrast, the maximum 
reversed velocity in the Gault NACA-0010 airfoil case, where the impact of 
windward differencing is greater, is u/u e = -.28. 

The predicted iso-velocity contours using the windward differencing 
scheme are compared in Fig. 13 with the experimental contours measured by 
Gaster. The iso-velocity contours predicted using the FLARE approximation 
are nearly identical to those shown for the windward scheme. Although good 
agreement was obtained with the experimental pressure distribution and the 
separation and reattachment locations as shown in Ref. 21, the predicted 
flow field away from the wall differs substantially from the measured 
iso-velocity contours. It is observed in Fig. 13 that the analysis 
significantly underpredicts the overall boundary layer growth throughout 
the bubble. This difference is probably strongly influenced by the 
assumption of instantaneous transition at s/L = 1 .0 which causes a slope 
discontinuity in the predicted iso-velocity contours at this location. 
Clearly more work is needed to include a natural transition model in this 
analysis, which is addressed in the next section of this report. 
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Natural Transition/Turbulence Model 

Upon incorporation of the McDonald-Fish-Kreskovsky (MFK) model into 
the ALESEP code, a series of direct boundary layer calculations using the 
present finite difference boundary layer procedure were performed to 
demonstrate the capabilities of the model for attached flows. For 
separated flow, the modified McDonald-Fish (MF) model was used in the 
identical manner as Briley and McDonald (Ref. 15) did in their study of 
transitional separated flow on airfoils. In all of these calculations a 
total of 100 grid points were placed across the boundary layer with the 
minimum grid spacing located at the wall. 

Attached Flow 


Figure 14 shows the predicted skin friction and displacement thickness 
with the experimental data taken by Blair (Ref. 37) for transitional, zero 
pressure gradient flow over a flat plate at a high free stream turbulence 
level. The free stream velocity was 100.0 ft/sec and the unit Reynolds 
number was 5.105 x 10^. A total of 54 grid points were distributed 
unevenly in the streamwise direction with the minimum spacing located at 
the leading edge of the flat plate. In this case, the free stream 
turbulence decayed from 7 percent at the initial station to 3.6 percent at 
x = 94 inches. Figure 14 shows that the calculated transition region is 
between x = 2.2 and x = 6.2 inches. Since no detailed experimental data 
exists in the transition region of the flow due to the thinness of the 
boundary layer, this case demonstrates the accuracy of the boundary layer 
scheme in turbulent flow using the MFK turbulence model. This case is 
typical of the high free stream turbulence levels encountered in 
turbomachinery flows. In external aerodynamics, which is the focus of the 
present investigation, the free stream turbulence levels are generally 
substantially lower. 

By removing the turbulence generating grids upstream of the flat plate 
test section, Blair (Ref. 37) measured the boundary layer characteristics 
for a free stream turbulence level of 0.25 percent. The free stream 
velocity and Reynolds number were held the same as in the previous case. 

The computational grid defined in the previous example was used in the 
predictions for this case as well. Figure 15 shows the comparison between 
the results obtained with the present direct finite difference boundary 
layer solution technique using the MFK turbulence model and the experimen- 
tal data. It is observed that in this case, the predicted location of 
transition is significantly further downstream than that observed 
experimentally. Transition is predicted to occur between x = 37.2 and 55.2 
inches whereas transition was measured to occur between x = 27.4 and 36.0 
inches. Delayed transition was also predicted by the ABLE (Analysis of 
Boundary Layer Equations) finite difference direct boundary layer solver 
(Ref. 38) using the McDonald-Fish-Kreskovsky transition/turbulence model 
for this low free stream turbulence case. The accuracy of the present 
results in the laminar and turbulent regions is demonstrated in Fig. 15 
w ^th the inclusion of the known analytical behavior in these respective 
regions for the skin friction and displacement thickness for a flat plate. 


20 



Excellent agreement is observed in these two regions of the flow thereby 
confirming the basic approach for fully laminar and turbulent regions while 
highlighting the transition region shortfall. 

The calculations presented in Figs. 14 and 15 were for zero pressure 
gradient flow over a flat plate. By placing a contoured body opposite to a 
flat plate section, Sharma, Wagner, Edwards and Blair experimentally in- 
duced the boundary layer edge velocity distribution shown in Fig. 16(a) on 
the flat plate measuring surface. This experiment was conducted in the 
UTRC Boundary Layer Channel to investigate the boundary layer chracteris- 
tics in a flow with a pressure distribution which simulates that in the 
leading edge region of an airfoil. The inlet to exit velocity ratio was 
0.37 and the inlet free stream turbulence level was 7 percent. ^ The exit 
Reynolds number based upon a 36 inch test section was 6.0 x 10 . A com- 
parison between the calculated displacement thickness and skin friction of 
the direct finite difference boundary layer analysis and the experimental 
data is shown in Figs. 16(b) and 16(c). A total of 200 points were distri- 
buted evenly along the flat plate in the boundary layer calculation. Local 
turbulence levels, T y , were calculated from the inlet free stream turbu- 

lence and the local to upstream velocity ratio based upon the assumption of 
frozen turbulence. Although there is no experimental data in the transi- 
tional region, the agreement between prediction and experiment in the 
laminar and turbulent regions indicates that the predicted transition 
region is relatively close to that of the experiment. Transition was 
predicted to exist between x = 24.24 and 29.05 inches with the MFK model. 

A non-interacting inverse boundary layer calculation was run on this 
case to demonstrate the agreement between the direct and inverse boundary 
layer results using the MFK turbulence model. In order to be consistent 
between the direct and inverse calculations, the computed values of the 
perturbation mass flow parameter, m, resulting from the direct finite 
difference boundary layer calculation described above, were used as boun- 
dary conditions in the inverse calculation. The inverse boundary layer 
calculation reproduced the boundary layer edge velocity distribution shown 
in Fig. 16(a) demonstrating that the MFK turbulence model has been 
correctly implemented for inverse boundary layer calculations. 

This completes our calculations for attached boundary layer flows. 

These three cases demonstrate that the MFK transition/turbulence model 
gives good predicted results for flows with high free stream turbulence 
levels. However, for low levels of free stream turbulence, the MFK 
transition model predicts transition too far downstream resulting in a 
substantial difference between the predicted boundary layer quantities and 
the corresponding experimental measurements. This same weakness in the 
transition model is carried over to separated flows as will be demonstrated 
in the next section. 

Separated Flow 

Several calculations for transitional separation bubbles with low 
free stream turbulence levels have been attempted with the McDonald-Fish (MF) 
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turbulence model, modified for separated flow in the same manner as Briley 
and McDonald (Ref. 15) suggested. Among these calculations are the Gaster 
case, discussed previously, and the zero angle of attack flow over a NACA 

66 3 -018 airfoil experimentally studied by Gault (Ref. 3). This latter case 
was computed by Briley and McDonald. 

In the present calculations, it has been found that the onset of 
transition is predicted too far downstream with the MF transition model. 

The result has a catastrophic effect on the calculation of transitional 
separation bubbles since the predicted bubble size is now much too large 
with reattachment frequently not occurring thereby leading to massive 
stall. In contrast, Briley and McDonald successfully used the MF 
transition model in their time dependent Navier-Stokes /boundary layer 
analysis of the Gault mid— chord transitional bubble case with none of the 
difficulties which have been encountered in the present study with this 
transition model. The difference in the present calculations and those of 
Briley and McDonald is unresolved at the present time; there are signifi- 
cant differences in these two inviscid— viscous interaction formulations but 
whether or not these differences affect the calculation of the turbulent 
kinetic energy equation in the transition region will require further 
study. The delay in transition predicted in the flat plate flow with low 
free stream turbulence, which was discussed earlier, clearly points to the 
need for further work in adapting this transitional model to the low free 
stream turbulence flows in external aerodynamics. It has been found that 
reattachment can be forced to occur in the analysis of separated flows with 
the MF transition model by artificially increasing the free stream turbu- 
lence level, however. This is demonstrated using the flat plate experiment 
of Gaster. When the calculation was performed for the free stream turbu- 
lence level of 0.25 percent reported by Gaster, the ALESEP code with the MF 
model predicted that the flow separated from the flat plate with transition 
and thus reattachment not occurring in the region prior to the downstream 
end of the calculation located at s/L =1.5. By artificially increasing 
the free stream turbulence level to 4.0 percent, the flow was made to 
reattach within a reasonable length as shown in Fig. 17 where the predicted 
pressure, displacement thickness, and skin friction with this level of free 
stream turbulence are predicted. Local edge turbulence levels were calcu- 
lated using the local to upstream velocity ratio using the assumption of 
frozen turbulence. A further increase in turbulence level to 4.5 percent 
results in a significantly smaller bubble with the pressure distribution 
not showing the usual pressure plateau region characteristic of separated 
flows. To further demonstrate the sensitivity of these calculations to the 
prescribed free stream turbulence level, it was found that a free stream 
turbulence level of 3 .5 percent results in a bubble which has more than 
doubled in size compared to that at 4.0 percent. 

In order to demonstrate the accuracy of the turbulence modeling in the 
MF model, transition can be forced to occur by simply varying the struc- 
tural coefficient, a^ , from 0.0 in laminar flow to 0.15 for turbulent flow 
over a specified transition length. The variation between these levels is 
specified with the Dhawan and Narasimha intermittency distribution. The 
angle of attack for this case is 0.0 degrees and the chord Reynolds number 
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is 2.0 x 10^. The reference pressure distribution was taken to be the 
experimental distribution which was obtained at high Reynolds number (Re = 
1 x 10 7 ) in which transition naturally occurred before laminar separation 
could take place. A direct boundary layer calculation was performed from 
the leading edge of the airfoil to s/c = 0.5 using the reference pressure 
distribution as the edge boundary condition. The velocity profile was 
taken at the s/c =0.5 station from the direct calculation to be used as 
the initial profile for the interacting calculation. The reference 
displacement thickness was held constant at that value predicted by the 
direct boundary layer calculation at s/c = 0.5. A total of 99 points were 
distributed evenly between s/c = 0.5 and s/c = 0.99 in the interaction 
calculation. The local edge turbulence level was calculated from a 0.2 
percent free stream turbulence level and the local to upstream inviscid 
veloctiy ratio using the frozen turbulence assumption. Transition was 
prescribed to occur between s/c = .693 and s/c = .703. Figure 18 shows the 
results of the interaction prediction using this forced transition model on 
the NACA 662 -OI 8 airfoil tested experimentally by Gault (Ref. 3). The good 
agreement between the predicted pressure distribution and Gault* s 
experimental pressure data shown in Fig. 18(a) demonstrates that the 
McDonald-Fish turbulence model may be successfully used in the interaction 
calculation providing that transition is predicted correctly. Figures 
18(b) and (c) show the predicted displacement thickness and skin friction 
distributions for this case. The predicted boundary layer velocity 
profiles at six locations along the airfoil surface using both the FLARE 
and windward schemes are compared with experimental profiles taken by Gault 
in Fig. 18(d). For this case where the maximum reversed flow velocity was 
u/u e = -.08, the FLARE and windard schemes produced nearly identical 
results. The predicted profiles differ from the experimental profiles in 
the transition region but agree fairly well in the laminar and turbulent 
regions of flow. As previously shown by Carter and Vatsa (Ref. 21), the 
predicted results using the interaction technique is quite sensitive to the 
location of transition. It is shown in Fig. 19 that the sensitivity 
remains unchanged using the transition model described above. A slight 
change in the intermittency distribution shown in Fig. 19(a) leads to a 
major change in the predicted results as shown by a comparison of the 
experimental velocity profiles at s/c = .725 with the predicted profiles in 
Fig. 19(b) using these intermittency distributions. A reduction in transi- 
tion length of only 0.25 percent chord results in a 22 percent shorter 
separation bubble and thus an attached turbulent profile at this location. 
Figure 19 dramatically demonstrates the accuracy which is required of a 
transition model in order to accurately predict the flow in a transitional 
separation bubble. 


CONCLUDING REMARKS 

The development of an improved airfoil transitional separation bubble 
analysis has continued through the inclusion of a proper finite difference 
technique in the reversed flow region and the incorporation of a natural 
transition/turbulence model to predict the onset of transition within the 
separation bubble. This introduction of a windward finite difference 
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procedure into the present inviscid-viscous technique has been found to be 
stable and inexpensive in terms of additional computational time over that 
with the Reyhner-Flugge Lotz (FLARE) approximation. It has been found for 
flows with reversed flow velocities up to 28 percent of the boundary layer 
edge velocity that comparable results for the surface pressure, skin 
friction, and displacement thickness are predicted by the windward and 
FLARE schemes. However, the use of the more accurate windward difference 
scheme has revealed a major structural change in the interior of the 
separation bubble. A counter-rotating region under the primary transitional 
separation bubble can emerge with the use of windward differencing which 
has not previously been observed with the FLARE approximation. Numerical 
tests were performed to indicate that the existence of this secondary 
structure is not sensitive to grid spacing. The differences in bubble 
structure and boundary layer characteristics between the windard and FLARE 
schemes are proportional to the size of the reversed flow region and the 
magnitude of the reversed flow velocity within the bubble. Convergence 
studies also indicate that with both the windward and FLARE schemes, the 
present inviscid-viscous coupling procedure is capable of reducing resid- 
uals to desired levels. The McDonald-Fish-Kreskovsky natural transition 
turbulence model has been incorporated into the interaction code for the 
prediction of transition within the separation bubble. A number of cases 
are presented which demonstrate that this model predicts transition too far 
downstream in separated and low free stream turbulent flows. Sensitivity 
studies are reported which indicate that modifications to this natural 
transition model are required for these flows. 
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APPENDIX A 


Analysis of the Inviscid Flow Over an 
Airfoil with a Small Protrusion 

In this appendix the method of matched asymptotic expansions is used 
to analyze the inviscid flow over an airfoil on which a small protrusion 
occurs. Figure A-l shows a schematic diagram of the airfoil with the 
protrusion details enlarged so that the length scales which characterize 
this bump can be clearly shown. First, a general formulation of the 
overall inviscid flow will be presented followed by a detailed description 
of the outer flow over the airfoil and inner flow in the immediate vicinity 
of the protrusion. Matching of these two flows through second order will 
be discussed with the end result being a mathematical formulation of the 
lead order perturbed inviscid flow near the protrusion which is used in 
the main text in the local inviscid-viscous interaction analysis. The 
present analysis is conducted only for incompressible flow; however, the 
inclusion of compressibility effects should straightforwardly follow the 
formal approach taken here. 

The incompressible, irrotational flow over an airfoil can be described 
by the full potential equation 

(tT® £ ) £ + 0 (A-l) 


where £ and n are orthogonal curvilinear coordinates oriented to the smooth 
airfoil sur f ace as shown in Fig. A-l. The subscripts £ and n denote diffe- 
rentiation in the respective directions. For convenience the origin of the 
£ ,n-coordinate system has been chosen at the center of the protrusion as 
shown in Fig. A-l. The metric coefficient H is defined as H = 1 + where 
k ie the curvature of the airfoil surface (without the protrusion). The 
tangential and normal velocity components, U and V, respectively are 
related to the velocity potential $ by 


U 





The surface boundary condition is given by 



(A-2) 



(A-3) 


which is imposed at n o = 6T(£) and requires the flow to be tangent to the 
airfoil, including the protrusion. Equation (A-3) can also be expressed in 
terms of the velocity potential which results in 


$ 


7 ? " 




at 7 ) = i 7 o 


(A-4) 
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In the present analysis, it is assumed that the streamwise length of the 
protrusion is e where e is asymptotically small compared to the chord of 
the airfoil. Hence, it is seen that over most of the airfoil the usual 
Neumann condition, (£,0) = 0, is the required boundary condition. In 
the far field, the flow is required to return to uniform flow which is 
stated by 

<t>£ — h cos 9 I 

/ as £ , 17 00 (A- 5 ) 

4^ — - sin 6 l 

where 0 is the airfoil surface angle measured with respect to a reference 
axis . 

It is convenient to break the analysis into an outer or global description 
of the flow, and an inner or local description which is applicable in the 
immediate vicinity of the protrusion. The relative length scales of the 
bump are chosen also as to be representative of those of viscous displace- 
ment thickness induced by a transitional separation bubble. In our 
previous investigation (Ref. 21) it was observed that the height/length 
ratio of the displacement surface produced by a transitional separation 
bubble was always very small. Therefore, in the present study it is 
assumed that 


€ 




(A-6) 


and hence, the protrusion height, 6, is taken to be asymptotically small 
compared to its length, e, in the limit of a vanishing protrusion. We now 
proceed to establish the outer and inner region expansions in terms of the 
protrusion length, e. 


Outer Region 

The outer potential is expanded in terms of the gauge functions A n to 

give 


<t, ° (C , V ) = A, (i'*)) + ^2 ( € ) *** 2 (C'V) + (A-7 ) 

where each of the are assumed to be 0(1). Substitution of Eq. (A-7) 
into Eq. (A-l) and taking the limit as e ^ 0 it is seen that each of the $ 
satisfies the original partial differential equation. 

hr *k) t < A - 8) 

The curvature of the airfoil, k, is assumed to be 0(1), that is, as e 0 
the airfoil curvature remains finite. 
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We now substitute the expansion for the outer potential into the 
surface boundary condition given by Eq. (A— 4). Since we are assuming the 
protrusion is thin, that is, its height is small compared to its length, 
the boundary condition can be transferred via Taylor series from the 
protrusion surface n = ri Q to n = 0, the original airfoil surface. The 
surface boundary condition now becomes 


% (£,o) - -f- T' 


<t>°, le C\ A 2 *° 


>,£ U, o) - ^2. U ,o) 


+ 8 t| 4>° 


7J7J 


U,°) + <J>° 27?7? (C.o) 


A, 



U,0) 


(A- 9) 


+ -— tt ' \<t>° liv U,o) - (£, 0 )| + ... =0 

In Eq. (A-9), T' denotes dT/d£ which is an 0(1) quantity with £ defined as 
^/e. Equating 0(1) terms to zero yields 


tfr) (C.O) - 0 


(A-10) 


which is the usual surface boundary condition for the flow over an airfoil 
with no protrusion. 


Substitution of the outer expansion in Eq. (A-7) into the far field 
boundary condition yields 


. 0 , A 2 .0 

4> ,£ + —— <J> 


I 


2 C 


H COS 6 


<t>? 

1 V 


a 2 

a, 


sin 9 


as £, 


co 


(A-ll ) 


As e 0, it is concluded that A^ must remain finite and hence. A,, is set 
to unity for convenience. At this point, the description of the first 
order outer problem is complete and it is seen that it is simply the flow 
over the airfoil without a protrusion. 

Inner Region 

In the inner region, it is convenient to define stretched coordinates 


X = 



(A-12) 
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so that in these coordinates the protrusion length is 0(1) and the height 
is 0(5/e). Equation (A-6) shows that even in these coordinates the 
protrusion appears as a small disturbance to the local airfoil flow. 

The potential function for the inner region is expanded in terms of e as 


(£,7?<€) = 8, («) 4>[ (X,Y) + 8 2 U) <P Z (x,y) + 


(A-13) 


Substitution of this expansion into the governing equation (A-l) yields 


^IXX + ^lYY + ^2XX + ^2YY^ ~ e 


(A-14) 


+ € (^y</»; y ) y - 


8j € 


( K Y gx ^ x " ^ KY< £zY^Y + *** 


Equation (A-14) shows that in the limit of e -*■ 0 the first order inner 
potential is governed by Laplace's equation 


6 1 + <± ' = 0 

' I X X ' I YY U 


(A-l 5) 


The surface boundary condition, Eq. A-4, is expressed in terms of the inner 
region potential as 


= 

~ V 


_l_ 

H 


< 4>i 


(A-16) 


which is implemented at Y = Y q = <5/e T(X) where T(X) is an 0(1) function. 
Since the protrusion is treated as asymptotically thin, even in inner 
coordinates the surface boundary condition can be transferred to Y = 0 with 
a Taylor series expansion. Substitution of the inner expansion for the 
potential function Eq. (A-13), into the surface boundary condition, Eq. 
(A-16), and the use of a Taylor series yields the following relation to be 
satisfied at Y = 0: 


K (x .°> + -fr <#>ir < x '°> + -r I T <#• 1YV 

<X,°) - t'4>2X ( X ,0) I - TT' <£ 


(X , 0 ) - T>j x (X.O) 


hi 8 f 

8, ~ [ 


T 4> 


2 YY 


I XX 


(A- 17) 


+ 28 ^YT' <£ j x (X,0) + • • ♦ =0 

As e ■> 0 the 0(1) surface boundary condition is 


<*>/ Y (X,0)rO 


(A-18) 
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and thus we see that to first order neither the inner or outer solution is 
affected by the presence of the protrusion for the limit considered here 
(6/e « 1). 

Matching the 1-term inner limit of the outer expansion with the 1-term 
outer limit of the inner expansion yields 

€X<D° 1X (0,0) = 8, <Pj (00, CD) (a _ 19) 

In order for these solutions to match we set 5. = e. The solution for the 
inner first order potential function is now established by Eqs. A-15, A-18, 
and A-19 and is simply 


4>;(X,Y) = U°(0,0)X (a _ 2i 

where U°(0,0) = ^^(0,0) is the airfoil surface speed evaluated at the 
protusion location, 5=0, n = 0, on the airfoil surface. So we see that 
the influence of the protrusion is relegated to the second order problem 
which will be discussed next. 

Second Order Analysis 


The governing equation for the second order outer problem has already 
been given as Eq. (A-8) with n = 2. The surface boundary condition, Eq. 
(A-9) is written as 


(£. 0 ). 

iT ) 


€ &2 


V <t> 


• £ 


(£.o) 


A2 


T <x> c 


7J7J 


(£,o) 


= o 


(A-21 ) 


This relation shows that the inner boundary condition on the second order 
outer potential is determined by the choice of the gauge function, If 

we choose A 2 - 5/e, the height to length ratio of the protrusion, then the 
inner boundary condition for the second order outer problem will be the 
usual small disturbance boundary condition given by the first two terms in 
Eq, (A-21). However, since we want this boundary condition to be met by 
the inner solution, A 2 is chosen such that that is 


thereby resulting in 



(A-22) 


<x>2^(£,o) =o 


(A- 23) 


Since the outer boundary condition for the second order outer potential is 
homogeneous as deduced from Eq. (A— 11), then the solution for the second 
order outer potential which satisfies Eq. (A- 23) and the governing equation 
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is the trivial solution = 0. This result is expected since the 

first order inner solution was found to be unaffected by the presence of 
the protrusion. 

The second order inner problem is now examined and will show for the 
first time the explicit appearance of the protrusion characteristics. The 
governing equation, (Eq. (A-14), is rewritten as follows with the first 
order inner solution inserted: 


+ + •••=0 < a - 24 > 

The gauge function, 6 £> i s chosen so that 62 > e 3 which results in 

<#>Ux + =0 (a ' 25) 

The surface boundary condition, Eq. (A-17) is rewritten as 

<^Y (X * 0) " r (x,0)+ -J- | <^)^ YY (X,°) - T 7 <£j, x (x,0)j 


2 Sc 
+ S 2 


/cYT 7 




+ 


= 0 


(A- 26 ) 


Since we want the inner solution to describe the flow over the protrusion 
the gauge function, 6 £ > is set equal to the protrusion height, 6 , which 
results in 


<^2y (x .°) = U° (0,0) T' 


(A-27 ) 


Equation (A-27) is the usual small disturbance surface tangency condition 
which in this case is that of a thin protrusion subjected to the local 
velocity of the outer solution evaluated at the protrusion origin. 

The outer boundary condition on the second order inner potential is 
determined by matching the inner and outer expansions. The matching 
condition between the 2 -term inner and 2 -term outer expansions leads to 

44 (00 * CD ) = T5" [ X 2 $?££ (0,0) + XYG 0 ,^ (0,0) + Y 2 4> C i 7? ^(0,0)J(A-28) 


The least restrictive condition that can be imposed at this point is to 
choose the limit 


lim 


€ — 0 


8 


= o 


(A-29) 
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which results in 


$2 (® .00) =0 (A-30) 

At this point two-term inner and outer expansions have been derived 
which describe the local disturbance due to a protrusion placed on the 
airfoil. For the particular limit where the bump height, 6> relative to 
the length, e, satisfies the inequality < 6 < e, it has been determined 
from matched asymptotic expansions that to lowest order the disturbed flow 
over the protrusion is governed by Laplace* s equation with the usual small 
disturbance boundary condition at the airfoil surface and uniform flow at a 
large distance from the protrusion. This uniform flow velocity is that of 
the first order outer solution over the airfoil evaluated at the protrusion 
origin. It is observed that to second order, the inner solution is 
independent of the airfoil curvature. Further examination shows that ex- 
plicit dependence on the airfoil curvature does occur in the third order 
inner solution. It is deduced from this analysis that to second order the 
outer solution is unaffected by the presence of the protrusion. 

It is of interest to compare the particular protrusion height to 
length ratio chosen in the present analysis with that inherent in triple 
deck theory for strongly interacting laminar flows. Triple deck theory is 
a multi-layer asymptotic analysis developed by Stewartson and Williams 
(Ref. 39) in which the wall layer, of length, = Re~3/8 and height, = 
Re~ ' , plays a role analogous to that of the bump in the present study. 
Triple deck theory has been successfully used in numerous strongly 
interacting viscous flow applications to delineate the mathematical 
structure as the Reynolds number tends to infinity. Even though triple 
deck theory has been limited thus far to laminar flows and the present flow 
under study is transitional, nonetheless it is of interest to determine as 
to whether or not triple deck theory overlaps with the present analysis. 

In the present work we are assuming that 6 > w hich is equivalent to 
stating that the bump height to length ratio is 


> length 

len S th bump 


(A-31 ) 


In triple deck theory the inner deck dimensions yield 


height 

length 


triple 

deck 


= length 


2/3 


(A-32) 


Since the length of the inner deck is vanishingly small as Re -> «>, it is 
seen that the triple deck height to length ratio satisfies Eq. (A-31) and 
hence, we can state that the particular protrusion dimensions chosen in the 
present analysis are inclusive of those contained in triple deck theory. 
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LEADING EDGE SEPARATION BUBBLE 
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Fig. 1 Schematic diagram of airfoil laminar-transitional separation bubble and pressure distribution 
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Fig. 2 Local interaction region coordinate system 
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Fig. 3 Finite difference cell structure for the boundary layer equations — attached flow. 


37 


4 i.j CONTINUITY (STREAM FUNCTION) 

^ i.j-1 




i-1 J 


e i.j + i 

• MOMENTUM AND ENERGY 

'J i + 1,j 
© i.j-i 


£ 


Fig. 4 Finite difference cell structure for boundary layer equations — reversed flow. 
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Fig. 5 Inviscid-viscous iteration procedure 
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Fig. 6 Comparison of results for windward and FLARE differencing — NACA-0010 
airfoil (modified). 

(a) Pressure distribution. 
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Fig. 6 Comparison of results for windward and FLARE differrencing — NACA-0010 
airfoil (modified). 

(b) Skin friction 
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Fig. 7 Transitional separation bubble streamline pattern — NACA-0010 
airfoil (modified). 

(a) Windward differencing 
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Fig. 8 Velocity profiles from windward and FLARE differencing — NACA-0010 airfoil (modified). 



Fig. 9 Effect of streamwise mesh size on skin friction from windward differencing. 
NACA-0010 airfoil (modified) 
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Fig. 10 Global convergence history — NACA-0010 airfoil (modified). 
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Fig. 11 Comparison of windward and FLARE differencing for the Gaster experiment, 
(a) Skin friction. 
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Fig. 12 Comparison of velocity profile at s/L = 1.0125 for windward and FLARE 
differencing for Gaster case. 
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Fig. 13 Iso-velocity contours for Gaster case. 





0.010 


0.008 


u M =100 ft/sec 
Re^/in =5.105 x10 4 

Tu INLET = 007 


PRESENT ANALYSIS 

□ EXPERIMENT (BLAIR) 


Cn 

N3 



0.002 


0 


1 


0 


20 


40 


60 


80 


100 


s (in.) 


Fig. 14 Flow over flat plate at high free stream turbulence level 
(a) Skin friction. 
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Fig. 14 Flow over flat plate at high free stream turbulence level, 
(b) Displacement thickness 
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Fig. 15 Flow over flat plate at low free stream turbulence level, 
(a) Skin friction. 




s (in.) 


Fig. 15 Flow over flat plate at low free stream turbulence level 
(b) Displacement thickness 




Fig. 16 Flow over flat plate with pressure 
(a) Boundary layer edge velocity 




Fig. 16 Flow over flat plate with pressure gradient, high free stream turbulence level 
(b) Displacement thickness 
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Fig. 16 Flow over flat plate with pressure gradient, high free stream turbulence level 
(c) Skin friction 
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Analysis of Gaster case with McDonald-Fish natural transition model, 
(b) Displacement thickness. 
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Fig. 17 Analysis of Gaster case with McDonald-Fish natural transition model, 
(c) Skin friction 
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Fig. 18 Predicted results for NACA 663-018 airfoil 
(b) Displacement thickness. 
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Fig. 18 Predicted results for NACA 663-018 airfoil, 
(c) Skin friction. 



Fig. 18 Predicted results 
(d) Velocity profil 
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Fig. 19 Sensitivity of computed results to transition location on NACA 663-018 airfoil 
(b) Velocity profiles at s/c =0.725 





Fig. A-1 Airfoil with a Small Protrusion 
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